• We are often interested in subdividing or segmenting the brain into meaningful biological regions of interest (ROIs) for an analysis.
• Examples: tissue segmentation, segmentation of gray matter structures, segmentation of pathology (MS lesions, tumors, …)
• We will perform 3-class tissue segmentation in R using fslr and ANTsR:
  • Cerebrospinal fluid (CSF)
  • Gray Matter (GM)
  • White Matter (WM)

• Let’s read in the training T1 and brain mask for subject 05 (not 01!).

library(ms.lesion)
library(neurobase)
all_files = get_image_filenames_list_by_subject(
  group = "training", 
  type = "coregistered")
files = all_files$training05 # NOT training subject 1!
t1 = readnii(files["T1"])
rt1 = robust_window(t1)
mask = readnii(files["Brain_Mask"])

run_mask = t1 > 100
dd_orig = drop_empty_dim(run_mask, keep_ind = TRUE)

• Many tissue class segmentations are based on k-means clustering.
• These methods can be skewed by large outliers.

hist(t1, mask = mask, breaks = 2000); text(x = 600, y = 3000, '"outliers"')

ortho2(rt1, t1 > 450, xyz = xyz(t1 > 450)) # xyz - cog of a region

t1[ t1 < 0 ] = 0
t1 = mask_img(t1, mask)
rt1 = robust_window(t1)

hist(rt1, mask = mask, breaks = 2000); 

• FAST is based on a hidden Markov random field model and an Expectation-Maximization algorithm (Zhang, Brady, and Smith 2001).
• It jointly produces a bias field corrected image and a probabilistic tissue segmentation.
• More robust to noise and outliers than finite mixture model-based methods that do not incorporate spatial information.

The fslr function fast calls fast from FSL. The --nobias option tells FSL to not perform inhomogeneity correction (N4 already performed in ANTsR).

t1file = files["T1"] 
t1fast = fast(t1, 
              outfile = paste0(nii.stub(t1file), "_FAST"), 
              opts = "--nobias")

FAST assumes three tissue classes and produces an image with the three labels, ordered by increasing within-class mean intensities. In a T1 image, this results in:

• Level 1: CSF
• Level 2: Gray Matter
• Level 3: White Matter

ortho2(rt1, t1fast == 3, col.y = alpha("red", 0.5), text = "White Matter")

ortho2(rt1, t1fast == 2, col.y = alpha("red", 0.5), text = "Gray Matter")

ortho2(rt1, t1fast == 1, col.y = alpha("red", 0.5), text = "CSF")

robust_fast = fast(rt1, # the robust_window(t1)
              outfile = paste0(nii.stub(t1file), "_FAST"), 
              opts = "--nobias")

• Overall the results look good
  • Not much difference after dampening outliers using robust_window

• Uses Atropos (Avants et al. 2011)
  • 3D K-means clustering + a Markov random field
• The extrantsr::otropos function works with nifti objects
  • calls ANTsR::atropos function

t1_otropos = otropos(a = t1, x = mask) # using original data
t1seg = t1_otropos$segmentation

ortho2(rt1, t1seg == 3, col.y = alpha("red", 0.5), text = "White Matter")

ortho2(rt1, t1seg == 2, col.y = alpha("red", 0.5), text = "Gray Matter")

ortho2(rt1, t1seg == 1, col.y = alpha("red", 0.5), text = "CSF")

• Overall the results do not look good
  • The k-means clustering is affected by large outliers
• We will try using robust_window

robust_t1_otropos = otropos(a = rt1, x = mask) # using robust
robust_t1seg = robust_t1_otropos$segmentation

double_ortho(rt1, robust_t1seg)

ortho2(rt1, robust_t1seg == 3, col.y = alpha("red", 0.5), text = "White Matter")

ortho2(rt1, robust_t1seg == 2, col.y = alpha("red", 0.5), text = "Gray Matter")

ortho2(rt1, robust_t1seg == 1, col.y = alpha("red", 0.5), text = "CSF")

• Overall the results look like they reasonably separate the classes
  • No ground truth
• Winsorizing large outliers aided the k-means clustering
  • Results much better than running Atropos on the raw data

We can create a table which will count the number of voxels in each category:

tab_fsl = table(robust_fast[ robust_fast != 0])
tab_ants = table(robust_t1seg[ robust_t1seg != 0])
names(tab_fsl) = names(tab_ants) = c("CSF", "GM", "WM")
tab_fsl

    CSF      GM      WM 
1551816 2632937 2932684 

tab_ants

    CSF      GM      WM 
1288931 2545481 3304745 

out_mask = robust_fast != 0 | robust_t1seg != 0 
table(robust_fast[ out_mask ], robust_t1seg[ out_mask ])

   
          1       2       3
  0   21720       0       0
  1 1267209  284607       0
  2       2 2257284  375651
  3       0    3590 2929094

By multiplying by the voxel resolution (in cubic centimeters) using the voxres function, we can get volumes

vres = oro.nifti::voxres(t1, units = "cm")
print(vres)

[1] 0.0001757813

vol_fsl = tab_fsl * vres
vol_fsl

     CSF       GM       WM 
272.7802 462.8210 515.5109 

vol_ants = tab_ants * vres
vol_ants

     CSF       GM       WM 
226.5699 447.4478 580.9122 

http://johnmuschelli.com/imaging_in_r